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SUMMARY 

\ 

The  model  described  in  this  paper  is  based  upon  numerical  integration  of  the  flow 
equations  which  simulate  the  water  movements  caused  by  tides,  wind  and  pressure  dif- 
ferences and  upon  numerical  integration  of  the  advective  diffusion  equation 
representing  Che  movements  of  salt  and  temperature.  The  model  neglects  vertical 
accelerations,  but  not  vertical  velocities,  and  includes  an  approximation  for  the  sub- 
gridscale  effects  by  introduction  of  mass  and  momentum  exchange  coefficients.  These 
coefficients  are  dependent  upon  the  Richardson  number  and  vertical  velocity  gradients. 
The  computation  method  is  tested  on  an  estuary  and  the  results  of  the  simulations  are 
presented  in  graph  form.  
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NOMENCLATURE 


■\x,\  - horizontal  momentum  dispersion  coefficients 
fc  = bottom  layer 
C = Chezy's  coefficient 
C*  = wind  stress  coefficient 
Dx,t>v  = horizontal  nass-heat  diffusion  coefficients 
! « Coriolis  parameter 
g = grav i t at iona 1 acceleration 
mi  (]  fk)  * layer  thickness 

i,i  = horizontal  grid  index 
k - vertical  grid  index 
1 = a layer  index 

m = a coefficient  related  to  the  vertical  mass  and  heat  diffusion  coefficient 
n = tine  level 
p = pressure 

r = a coefficient  related  to  the  vertical  momentum  diffusion  coefficient 
s = salinity 
T = temperature 

u,v,w  = respective  components  of  velocity 
wa  = wind  speed 

x,y,z  = Cartesian  coordinates,  positive  eastward,  northward  and  upward,  respectively 

■ « difference  operation 

' = water  surface  elevation  with  respect  to  the  reference  level,  e.g.,  mean  sea 
level 

■ = vertical  mass  diffusion  coefficient 

= vertical  thermal  diffusion  coefficient 
= a parameter  related  to  the  equation  of  state 

= vertical  momentum  diffusion  coefficient  under  nonisotropic  vertical  density 
grad ient 

= vertical  momentum  diffusion  coefficient  under  neutrally  stable  density 
gradient  (homogeneous) 

C = density 

~ = reference  density,  a constant 
. ' « departure  from  ~ depending  on  salinity  and  temperature 
' = shear  stresses 
; = wind  direction 
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INTRODUCTION 


In  many  estuaries,  three-dimensional  flows  occur.  The  fresh  river  water  often 
mixes  only  partially  with  denser  sea  water,  and  the  salinity  can  change  considerably 
over  the  depth.  The  currents  Induced  by  the  density  differences  are  superimposed  upon 
the  currents  caused  by  tides  and  wind,  causing  very  complicated  current  patterns. 

The  model  described  in  this  paper  can  effectively  simulate  these  flow  conditions 
by  making  the  assumption  that  the  vertical  accelerations  of  the  fluid  are  small.  This 
assumption  is  generally  justified.  Neglecting  the  vertical  accelerations  but  not  the 
vertical  velocities  makes  a mathematical  formulation  possible  which  permits  a direct 
numerical  integration  of  the  equations  representing  the  fluid  flow. 

The  continuous  fluid  properties  during  simulation  are  expressed  at  discrete 
points,  by  which  approach  certain  approximations  are  made;  the  approximations, 
however,  are  chosen  in  a manner  that  makes  certain  basic  physical  laws  of  the  fluid 
still  applicable.  For  example,  no  water  disappears  from  the  water  body  by  the 
numerical  procedures;  also,  energy  in  the  horizontal  motions  is  conserved,  disregarding 
the  losses  by  the  representations  for  turbulent  energy  dissipation  and  the  inputs  by 
wind . 


DYNAMIC  EQUATIONS  F 'R  TURBULENT  ESTUARINE  FLOW 

The  equations  of  horizontal  motion  for  an  incompressible  internally  source-free 
fluid  on  a rotating  earth  in  Cartesian  coordinates  with  the  z axes  positive  upward  are 
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If  the  vertical  acceleration  and  advection  are  neglected,  the  equation  of  motion 
in  the  vertical  direction  becomes  the  hydrostatic  equation 
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The  equation  of  continuity  is 
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The  equations  of  salt  and  heat  balance  are 
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and  the  equation  of  state 


p = p + p'  (s ,T) 


(3) 
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The  representation  of  the  estuarine  system  by  Eqs.  (1-7)  includes  no  internal 
sources.  In  the  equation  for  heat  balance  (Eq.  (6))  we  used  the  same  diffusion 
coefficient  as  in  the  equation  for  the  salt  balance  (Eq.  (5)),  which  means  that  we 
neglected  the  heat  fluxes  by  conduction  and  that  only  the  fluxes  by  the  turbulent  fluid 
are  considered. 
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The  model  described  In  this  paper  uses  a grid  system  with  equidistant  points  In 
the  horizontal  direction.  In  the  vertical  direction,  however,  the  possibility  of 
having  an  unequal  grid  distance  was  Introduced,  since  vertical  gradients  which  often 
occur  In  the  estuaries  can  then  be  represented  in  detail  without  using  a large  number 
of  grid  points  In  the  vertical. 

In  formulating  the  model,  we  introduce  the  concept  of  the  layer.  This  is  a 
geometrical  concept  only,  and  should  not  be  confused  with  layers  in  stratified  flow. 
The  layers  are  determined  by  the  introduction  of  a set  of  horizontal  planes  inter- 
secting the  body  of  the  fluid.  The  thicknesses  of  the  different  layers  are  not 
necessarily  the  same  for  all  layers. 

The  top  layer,  bounded  by  the  free  surface,  has  a time-variable  height.  In  the 
formulation  of  the  model  the  other  layers  may  be  intersected  by  the  bottom  and  have  a 
height  which  is  dependent  on  the  bathymetry.  Figure  1 illustrates  the  grid  system  and 
the  layers  of  the  model.  However,  in  our  experiments  we  did  not  vary  the  bottom 
layer. 


fig-  l--Locat ion  of  variables  on  the  vertical  grid 


In  the  computational  model,  the  origin  of  the  coordinates  is  taken  at  the  mean 
sea  level,  whereas  the  water  surface  z » C(x,y,t)  is  the  upper  boundary  of  the  system. 

The  description  of  the  finite-difference  equations  from  the  differential 
equations  was  accomplished  in  two  steps:  first,  the  equations  for  the  layers  were 
derived  by  vertically  integrating  the  variables  over  the  layer  thickness,  and  sub- 
sequently finite-difference  approximations  for  the  layer  equations  were  developed. 

We  have  shown  in  Refs.  1 and  2 that  the  vertically  integrated  momentum  equations  can 
be  written  for  each  layer  (Fig.  1)  using  u and  v for  layer  average  values  of  the 
velocity  components. 
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The  formulation  of  these  equations  of  motion  differs  from  the  more  classical 
representations  presented,  for  example,  by  Bowden  (Ref.  3)  and  Pritchard  (Ref.  4),  in 
that  our  equations  are  conceptionally  well  traceable  and  contain  the  derivatives  of 
products  of  variables  rather  than  products  of  variables  and  derivatives.  Except  for 
the  density  and  the  spatial  and  time  derivative,  no  dividers  are  present,  i.e.,  when 
finite-difference  approximations  of  these  equations  are  used,  the  number  of  divisions 
is  limited. 

The  heat  and  salt  balance  equations  can  be  presented  in  the  same  form  for  each 
layer  k: 
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From  the  equation  of  continuity  we  derived  for  the  time  derivative  of  the  water 
surface 
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and  for  the  vertical  velocity  at  the  interface  between  layer  k and  k - 1 
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The  horizontal  pressure  gradient  in  the  top  layer  was  approximated  by 
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where  p|  = layer  average  density  of  the  top  layer.  The  pressure  gradient  for  the  other 
layers  is 
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FINN  b-D I K FERENC K F0RMU1.ATI0N  AND  APPROXIMATION 


A space-staggered  grid  (Fig.  2)  was  selected  for  the  finite-difference  approxi- 
mation of  F.qs.  (8)  through  (17).  The  pressure  p,  salinity  s,  and  temperature  T are 
computed  at  integer  values  of  i,  j,  and  k.  The  velocity  components  u are  computed  at 
half-integer  values  of  i and  integer  values  of  j and  k,  and  the  velocity  components  v 
are  computed  at  integer  values  of  i and  k and  half-integer  values  of  j.  The  vertical 
velocity  components  w are  computed  at  integer  values  of  i and  j,  and  half-integer 
values  of  k (see  Fig.  3).  The  water  surface  elevation  r,  is  expressed  only  on  a two- 
dimensional  horizontal  grid  and  has  integer  values  of  i and  j. 


•-1  c 4-  S- 

i-1  1 -i  i 1*S  <*1 

Fig.  2 — The  location  of  u (— ) , v (i) 
and  other  parameters  (-*-)  in 
the  space-staggered  grid 


' i ' — 


, I 


(* 

—4 

^.5  ,0  ,T 

V 

i 

y 

*“  • _ 

~ t 

T 

lii Li  _ 

* 

h.  * l^4Ct) 


Fig.  3 — Relative  position  of  the 
variables  in  the  model 


A compact  notation  was  introduced  which  not  only  describes  the  finite-difference 
approximation  accurately,  but  also  gives  this  approximation  the  same  appearance  as  the 
differential  equations.  The  position  of  the  variables  in  the  horizontal  grid  is 
indicated  by  indices  1 , j , indicating  a position  iAx  and  jAy  from  the  origin  of  the 
coordinates,  with  i , j = 0,  ±4,  tl,  ±3/2. 

The  vertical  position  is  determined  as  to  its  location  in  the  center  of  the 
layer  numbered  from  the  top  with  integer  k = 1 , 2 , 3 . . . or  at  the  horizontal  inter- 
faces with  half-integer  values  k = 1/2,  3/2,  5/2  ... 

Time  is  determined  by  the  number  of  time  steps  (nAt)  from  the  reference  time,  with 
n an  integer  value. 

The  following  sum  and  difference  notations  are  adopted: 

F*  = '-3  j F [ f i + 4)  Ax , J Ay  ,kAz,nAt  ] + F[(i  - 4)  Ax , j Ay , kAz  ,nAt  ] J (18) 

5XF  = — j F[ (i  + 4)  Ax, jAy, kAz, nAt]  - F[ ( i - 4)  Ax , j Ay , kAz, nAt ] j (19) 


These  are  shown  only  for  x,  but  are  also  used  for  y,  z,  and  t. 
The  special  notation  used  for  shifted  time  levels  is 


F[iAx ,j Ay ,kAz, (n  + l)At] 

(20) 

F[iAx, jAy ,kAz, (n  - l)At] 

(21) 
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With  th-'st-  spu.e-timf  sum-differencing  operators,  the  finite-difference  approximations 
of  i£qs.  (8)  through  (17),  for  the  level  k with  layer-average  values,  take  the  following 
form : 


<5 1 ' ' = " ) 5xOl*u)  + ty^V^!  3t  11  0 


(22) 


6^0) 


(h  u u ) - 5 (hrv  tr  ) - h 5 (u  w ) 
x y z 


rr-x--xy  1 rx.  1 

+ fh  v - — h 5 p - — 

— y v r — ' 


Xr  — X 

0 


( “)  * ( “j 

' 'k-Mi  'k-s 


—x 

P 


-x 


4 ! hxA  5 u!  + 6 , h*A  " 5 u| 
x|  x x |_  yl  y y |_ 


at  i + *i,  j , k,  n 


(23) 


5 (iP'v)  = - 5 (hXu  71)  - 5 (Ir* v ) - h^  5 ) 


rT^xy  1 1 

- fh-u  - — Ir  5 p — 

— y y ~y 

0 0 


Kz)  * hr) 

' 'k+4  ' 'k-»s 


4 , Ir  A 4 v + 5 h"  A 5 v 
x|  X X |_  y I y y I _ 


at  i,  j + k,  n (24) 


t (hs)  = - ix(hAus")  - 4^(hJ'vs-v)  - h5z(wsA) 


+ 5 !h*D  5 s!  +5  Jh^D  5 - (<5  s J + lx 5 s ) 

xl  x x y|  y y |_  \ z \ z -/k_,. 


' k+4 

at  i,  j , k,  n 


(25) 


T^ThTT  = - 5x(h*uTX)  - 5y  (hVf6 7)  - h4z(wTZ) 


+ 5 ]h*l)J  tJ  + 5 |hyD. 4 T,(  - (k's  T \ + (k'S  T I 

xl  hx|_  yl  h v |_  \ z -/  , \ z -/ 


k-Hs  ' k-=s 

at  i,  j , k,  n 


(26) 


The  density  is  computed  using  s(i,j,k,n)  and  T(i,j,k,n)  according  to  the  equation  of 
state : 

o = [5890  + 38T  - 0.375T2  + 3s ] / [ (1 779. 5 + 11.25T  - 0.0745T2) 

- (3.8  + 0. 01T)s  + 0.698(5890  + 38T  - 0.375T2  + 3s)]  at  i,  j,  k,  n + 1 


(27) 


A graphical  representation  of  the  density  as  a function  of  temperature  for  water  with 
different  salinity  is  shown  in  Fig.  4.  The  finite-difference  equation  used  to  com- 
pute the  vertical  velocities  is 
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Fig.  4 — Equation  of  state  and  graphical  representations 
of  several  selected  ranges  of  s and  T values 

This  equation  is  used  for  the  bottom  layer  first  and  then  for  the  layer  above,  etc. 


The  horizontal  pressure  gradients  are  computed  first  in  the  top  layer: 

S P “ 8(0*6  C)  + *-S  h*5  P at  i + 4,  j , 1,  n + 1 (29 

A A X 
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Then  the  pressure  gradients  are  computed  with  increasing  k by  use  of 

«z(«xp>  " 86XPZ  at  i + *5,  j , k + p,  n + 1 (31 

<S _ <<5  p)  - g6  pz  at  i,  J + p,  k + >5,  n + 1 (32 
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Once  the  horizontal  pressure  gradients  are  known,  the  water  level  and  velocities 
can  be  computed  again  by  the  sequence  of  finite  difference  equations  (22)  through  (32) 

Special  procedures  are  required  for  the  computation  of  variables  on  or  near  the 
boundaries.  The  seaward  boundary  presents  particular  problems,  as  described  in  Ref.  2 
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VERTICAL  MASS  AND  MOMENTUM  EXCHANGE  COEFFICIENTS 

In  fluids  with  vertically  stable  or  unstable  density  gradients,  the  model  assesses 
the  vertical  momentum  exchange,  vertical  mass-heat  diffusion  coefficients  according  to 
the  following  Richardson  number  related  formulas  (Ref.  2).  For  the  vertical  momentum 
exchange : 
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For  the  vertical  diffusion-dispersion  ev  salt  and  heal: 
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where  m,r  » exponents,  equal  to  0.4  and  0.8,  respectively,  according  to  Mamayev  (Ref.  5). 

A quadratic  expression  was  used  for  surface  wind-induced  stress  and  bottom 
current-induced  drag. 
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No  thermal  or  salinity  fluxes  at  the  bottom  are  introduced,  and  the  exchanges  at 
the  bottom  can  be  set  to  zero. 

NUMERICAL  EXPERIMENTS 


To  illustrate  the  model's  operation,  we  present  two  sample  cases.  The  first  one 
involves  the  simulation  of  three-dimensional  wind— driven  circulation  in  a rectangular 
basin.  This  basin  was  3400  m long,  1400  m wide,  and  7 m deep.  This  rectangular  basin 
was  then  schematized  into  a grid  structure  of  17  « 7 * 7 with  Ax  » Ay  « 200  m and  Az  ■ 
1 m.  The  water  was  at  rest  at  the  beginning  (t  ■ 0).  A diagonal  wind  stress  of  10 
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dynes/cm^  (approximately  equivalent  to  40  mph)  was  applied,  and  the  simulation  was 
carried  out  to  2000  timesteps  with  At  « 5 sec.  Figure  5 shows  the  vertical  sections 
through  the  longest  axis  of  the  basin  In  a sequence  (a,b,c  and  d)  representing  the 
system’s  transient  response  under  wind  stress.  Each  graph,  with  the  simulation  time 
step  indicated  at  the  lower  left  corner,  shows  the  sequential  development  of  the  wind- 
induced  circulation  pattern  in  the  x-z  direction.  In  these  graphs  the  vertical  scale 
has  been  enlarged  200  times,  and  also  the  vertical  velocity  components  (i.e.,  w)  are 
multiplied  by  200  before  the  computation  of  the  x-z  velocity  vectors  at  each  grid 
point . 

In  the  beginning  only  a surface  current  is  set  up  by  the  wind  stress  (Fig.  5a). 
Figure  5b  shows  the  penetration  of  surface  current  to  the  deeper  layers  through  the 
viscosity  and  vertical  advectlon  effects,  and  vertical  velocity  components  at  two  ends 
of  the  basin.  Figures  5c  and  5d  show  the  generation  of  two  large  vertical  eddies, 
while  the  surface  current  gradually  penetrates  from  the  surface  layer  through  the 
bottom.  Figure  6a  shows  the  two  elongated  eddies  finally  breaking  up  and  merging  into 
one  large  basin-sized  eddy  forming  a layered  flow  structure.  At  the  75th  time  step 
(Figs.  6a  through  6h) , a quasi-steady  state  circulation  is  almost  established.  The 
effect  of  bottom  friction  becomes  more  pronounced  as  the  undercurrent  reaches  its 
maturity.  The  current  delimitation  plane  between  the  surface  current  and  the  return 
current  is  moved  upward  due  to  the  bottom  current's  retardation. 

Figure  7 shows  the  time  histories  of  the  multilayered  current  at  the  center 
section  of  the  basin  for  the  entire  simulation  period.  The  duration  of  the  simulation 
is  selected  approximately  equal  to  ten  times  the  basin  oscillation  period  of  the  basin 
(equal  to  approximately  164  time  steps).  In  the  first  few  hundred  steps  the  current 
field  becomes  established  and  the  basin-wide  circulation  cell  is  generated.  At  the 
beginning  of  the  computation,  wind  stress  suddenly  applies  energy  to  the  system.  Con- 
sequently, the  basin  also  starts  to  oscillate  and  the  velocities  belonging  to  this 
oscillation  are  superimposed  upon  those  of  the  wind-driven  circulation.  The  mean 
velocities  at  the  top  and  bottom  layer  are  shown  as  dashed  lines.  At  step  500  the 
wind  stress  is  removed  and  the  circulation  decays.  The  oscillation  of  the  basin,  how- 
ever, persists. 

The  Chesapeake  Bay  system  in  the  eastern  part  of  the  United  States  was  also 
selected  for  a numerical  experiment.  The  hydrodynamic  behavior  of  Chesapeake  Bay  is 
characterized  by  its  circulatory  responses  under  four  major  driving  and  modifying 
forces,  namely,  tide,  wind,  river  inflow,  and  Coriolis  forces.  The  geometry  of 
Chesapeake  Say  is  schematized  into  a computational  grid  network  of  40  x 100  x 6,  with 
Ax  = Ay  = 3000  m.  The  thickness  of  the  top  five  layers  is  4 m each,  with  the  last 
layer  being  8 m. 

Figure  8 is  a typical  simulation  output  map  showing  the  horizontal  circulatory 
pattern  for  the  top  layer  at  time  step  4000,  with  At  » 60  sec  (i.e.,  66.7  hr).  The 
ratio  between  the  length  of  a velocity  vector  and  a grid  size  for  this  particular  case 
is  30  cm/sec.  The  heavy  outline  denotes  the  water's  edge  at  the  lower  low  water, 
whereas  a thin  dash  line  delineates  the  20-ft-bottom  contour  line.  On  this  map,  the 
horizontal  velocity  field  is  represented  by  the  velocity  isollnes.  At  this  point  in 
time,  the  tide  level  at  the  open  boundary  is  near  its  low  slack  water  (see  input  tide 
at  the  lower  right  corner). 

In  the  upper  bay  water  flows  generally  to  the  north.  The  vertical  circulation 
pattern  at  a few  selected  sections  can  be  seen  from  the  vertical  current  profiles  in 
the  upper  right  corner  of  Fig.  8. 

CONCLUSIONS 

The  three-dimensional  flows  in  estuaries  with  nonlsotropic  density  distributions 
can  effectively  be  computed  according  to  the  finite  difference  method  described  in 
this  paper.  The  model  required  coefficients  used  in  describing  the  small-scale  mass 
and  momentum  exchanges.  These  coefficients  must  be  found  experimentally  from  field 
data. 
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Fig.  8 — A typical  current  map  showing  the  spatial  distribution  of  the  velocity 
field  at  2.0  m and  the  vertical  circulation  through  three  selected 
cross-sections  of  Chesapeake  Bay,  U.S.A. 
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